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The merger of two neutron stars launches a relativistic jet, which 
must be driven by a strong large-scale magnetic field. However, the 


magnetohydrodynamical mechanism required to build up this magnetic 
field remains uncertain. By performing an ab initio super-high-resolution 
neutrino-radiation magnetohydrodynamics merger simulation in full 
general relativity, we show that the «Q dynamo mechanism, driven by the 
magnetorotational instability, builds up the large-scale magnetic field inside 
the long-lived remnant of the binary neutron star merger. As a result, the 
magnetic field induces a Poynting-flux-dominated relativistic outflow with 
an isotropic equivalent luminosity of -10° erg s” and a magnetically driven 
post-merger mass ejection of -0.1 M,. Therefore, the magnetar hypothesis, 
in which an ultra-strongly magnetized neutron star drives a relativistic jet 
in binary neutron star mergers, is possible. Magnetars can be the engines 
of short, hard gamma-ray bursts, and they should be associated with very 
bright kilonovae, which current telescopes could observe. Therefore, this 
scenario is testable in future observations. 


The observation of GW170817/GRB 170817A/AT 2017gfo made binary 
neutron star mergers a leading player of multi-messenger astro- 
physics™. It revealed that at least a part of the origin of short, hard 
gamma-ray bursts and the heavy elements synthesized by rapid neutron 
capture (the r process) is binary neutron star mergers’ “°°. 

All the fundamental interactions play an essential role in binary 
neutron star mergers. Thus, the method chosen to explore them 
theoretically is a numerical relativity simulation that implements 
all the effects of the fundamental interactions. Numerical relativity 
simulations can qualitatively explain AT 2017¢fo, that is, the kilonova 
emissions associated with the radioactive decay of the r-process ele- 
ments” 8, However, a quantitative understanding of this event is still 
being developed. Moreover, there is no theoretical consensus about 
how the binary neutron star merger drove the short, hard gamma-ray 
burst? ”. 


A relativistic jet was launched from this binary neutron star 
merger, which was observed as a short and hard gamma-ray burst™*. 
Such relativistic jets are most likely driven by a magnetohydrodynam- 
ics process. This indicates that a binary neutron star merger remnant 
must build up a large-scale magnetic field with the dynamo to launch 
the jet”. However, the mechanism behind the large-scale dynamo in 
the merger remnant still needs to be clarified”. 

Also, the site for generating the large-scale magnetic field after the 
merger is a riddle. Is it inside a massive torus around a black hole that 
formed after a merger remnant collapsed into it? Or is it the long-lived 
remnant ofa massive neutron star”? Recent numerical simulations for 
a black hole and torus system, as a remnant of a binary neutron star 
merger or a black hole-neutron star merger, suggest that a large-scale 
magnetic field is established after a long-term evolution of O(1-10) s 
(refs. 20,27-29). Consequently, a relativistic jet is launched due tothe 
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Fig. 1| Overview of the magnetic field evolution in the remnant of the binary 
neutron star. Top, Electromagnetic energy as a function of the post-merger time 
for the total (solid), the poloidal (dashed) and the toroidal (dotted) components. 
The purple curves are for the simulation with AX ines: = 200 m. The inset shows 
how the amplification of the magnetic field due to the Kelvin-Helmholtz 
instability depends on the initial magnetic field strength and grid resolution. B14 
and B15 mean Bo max= 10” G and 10" G, respectively. Bottom, Magnetic field lines 
for a density of p <10” g cm? at t- t merger = 130 ms. The core of the hypermassive 
neutron star is shown for a density of p > 10° g cm”. 


Blandford-Znajek mechanism”. However, the long-lived remnant 
of a massive neutron star is more computationally challenging than 
the black hole plus torus system because the requirement to resolve 
relevant magnetohydrodynamical instabilities numerically is severe”. 
The physical mechanism for generating the large-scale magnetic field 
and launching the jet could differ in the two scenarios mentioned 
above. It could be observationally testable if we managed to build a 
long-lived binary neutron start model that generates a large-scale B 
field and launches a jet. 

We tackle this problem with a super-high-resolution neutrino- 
radiation magnetohydrodynamics simulation of a binary neutron star 
merger in full general relativity. 

We employ the latest version of our numerical relativity neutrino- 
radiation magnetohydrodynamics code”. We employ a static mesh 
refinement with 2:1 refinements to cover a wide dynamic range. For 
the simulations in this article, the grid structure consists of 16 Carte- 
sian domains. Each domain has 2N x 2N x N elements in the x, y and z 
directions, and we assume orbital plane symmetry. We set N = 361 and 
AX finest = 12.5 m. The grid resolution is the highest among the binary 
neutron star merger simulations”. The initial orbital separation is 
~44 km (Methods). 

We employ DD2 (ref. 33) as an equation of state for the neutron 
star matter and symmetric binary, which has a total mass of 2.7 Mo. 
With this set-up, a hypermassive neutron star that transiently formed 
after the merger will survive for >O(1) s (ref. 34). 

The purely poloidal magnetic field loop is embedded inside the 
neutron stars and has a maximum field strength of Bo,max = 105° G 
(ref. 35). It is much stronger than the 10’-10” G observed for binary 
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Fig. 2 | Generation of mean electromagnetic field. Mean and total 
electromagnetic energy in the region where the magnetorotational instability is 
active with p < 10"° g cm” as a function of the post-merger time. The blue, green, 
cyanand purple curves denote the mean poloidal, mean toroidal, total poloidal 
and total toroidal components, respectively, in the high-resolution simulation. 
The dashed curves are for the simulation with AX nese = 200 m. 


pulsars (ref. 36). The available computational resources limit us to a 
strong field strength and idealized topology. Nonetheless, it is natural 
to anticipate that in reality the amplification of the magnetic field leads 
toa high field strength (>10" G) in a short timescale after the merger 
(see below). 

Figure 1 plots the evolution of the electromagnetic energy as a 
function of the post-merger time. As reported in refs. 31,35, the elec- 
tromagnetic energy was exponentially amplified shortly after the 
merger due to the Kelvin-Helmholtz instability, which emerges at the 
contact interface when the binary merges. A part of the turbulent 
kinetic energy due to the instability is converted to electromagnetic 
energy”. Inthe inset, we generate the same plot for -1 < t ~ t merger $ 5 MS 
but different initial magnetic field strengths of Bo max = 10% G (red) 
and Bo max = 10 G (orange) while keeping the grid resolution. Also, 
we plot the results for the other grid resolutions Of AX fines: = 100 or 
200 m (brown and purple curves, respectively), while keeping 
Bo, max = 104° G. The figure clearly shows that the growth rate of the 
electromagnetic energy for 0 < t- tmerger £ 1 ms does not depend onthe 
initial field strength but on the grid resolution, as expected from the 
properties of the Kelvin-Helmholtz instability’”** (Extended Data 
Fig. 1). The post-merger amplification of the magnetic field due to this 
instability terminates when the shear layer is dissipated by the shock 
waves at t — tmerger * 2 MS. 

Att- Cnerger “ 5 ms, the electromagnetic energy temporarily settles 
to -3 x 10°° erg. However, the toroidal magnetic field (dotted curve) 
is subsequently amplified. In particular, its contribution to the total 
electromagnetic energy becomes prominent for t — t merger 2 20 ms and 
the growth rate is proportional to t°, which indicates that the mag- 
netic winding is efficient with a coherent poloidal magnetic field. The 
differential rotational energy of -1-2 x 10° erg of the remnant of the 
massive neutron star is the energy budget for the magnetic winding. 
Once the magnetic braking starts to work, the electromagnetic energy 
saturates around -5 x 10° erg. We also plot the result of the simulation 
with AXfinest = 200 m in Fig. 1. With this resolution, it is hard to resolve 
the Kelvin-Helmholtz and magnetorotational instability, particularly 
in the high-density region (Extended Data Figs. 1-3, and the section 
on convergence in Methods). Consequently, the amplification of the 
electromagnetic energy is less efficient in the low-resolution case. In 
particular, there is a striking difference in the poloidal magnetic field 
between the simulations with AX fines: = 12.5 and 200 m. Because the 
magnetic field that is amplified due to the Kelvin-Helmholtz instability 
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Fig. 3 |a dynamo inside the remnant of the binary neutron star merger. Butterfly diagrams at R = 30 km. Top left, Mean toroidal magnetic field By. Top right, Mean 
radial magnetic field Bg. Bottom left, Toroidal electromotive force Eg. Bottom right, Parameters a,, (orange) and correlation between Ep and By (blue). 


Table 1| The aQ dynamo period prediction and simulation 
data at several radii 


R(km) a,,(cms*) Q Shear k,(cm™) — Pmeory (S) Psim (S) 
(rads") rate 

20 -8.1x10° 4,025 q=1.00 6.3x10% 0.020 0.018 

30 -1.0x10 2,515 q=1.34 4.2x10° 0.021 0.018- 
0.024 

40 -1.0x107 1,688 q=1.44 3.3x10° 0.037 0.018- 
0.030 

50 -4.4x10° 1,200 q=1.50 2.6x10° 0.062 0.030- 
0.040 


Prheory, the aQ dynamo period; Psim, the Butterfly diagram period in the simulation. 


is randomly oriented”, there must be a mechanism to generate the 
coherent poloidal magnetic field. 

To make it clear which part of the merger remnant is responsi- 
ble for the generation of the coherent magnetic field, we evaluate 
the electromagnetic energy in the region where the magnetorota- 
tional instability is active, as defined by p < 10° g cm” and as shown 
in Fig. 2 (Extended Data Fig. 2). We also decompose the contributions 
from the mean poloidal and toroidal magnetic fields by taking an 
axisymmetric average. The electromagnetic energy in this region is 
subdominant compared to that contained deep inside the core region 
where p = 10° g cm” (see also Fig. 1). However, in the simulation with 
AXinest = 12.5 m, the mean poloidal field exponentially grows in the 
range 20 < t- tmerger £ 50 mS. At t — t merger * 20 ms, its contribution to 
the total poloidal field energy (solid cyan curve) is -1% and becomes 
~10% after the exponential growth. Such a clear exponential growth is 
invisible in the simulation with AX¢nes: = 200 m for t- tmerger £ 100 ms. 
Also, the mean poloidal field energy differs by one or two orders of mag- 
nitude in the simulations with AXfinest = 12.5 or 200 m. For the toroidal 
component, the mean and total field energy are of the same order of 
magnitude for t- tmerger 2 30 ms. We also confirm that the magnetoro- 
tational instability is responsible for generating the mean poloidal 
magnetic flux (Extended Data Figs. 4 and 5). These results indicate 
that the generation of a coherent magnetic field is triggered by the 
magnetorotational instability. 

The aQ dynamo, whichis driven by the magnetorotational instabil- 
ity” in the current context, is a potential mechanism for generating 
the large-scale magnetic field””. In the mean field dynamo theory, we 
assume that each physical quantity Qis composed of the mean field Q 


and the fluctuation q, that is, Q = Q + q. Thus, we cast the induction 
equation as 

0,B=Vx(UxB+é), (1) 
where € = u x b is the electromotive force due to the fluctuations, 
B = B + bis the magnetic field and U = U + uis the velocity field. Note 
that the magnetorotational instability-driven turbulence produces 
the fluctuations u and b. Inthe «Q dynamo, we express the electromo- 
tive force as a function of the mean magnetic field: 


Ēi = aB; + Bi(V x B), (2) 


where a, and £;are tensors that do not depend on B;. We calculate the 
mean field by taking the average in the azimuthal direction. 

In the presence of a cylindrical differential rotation, the simplest 
mean field dynamo is an «Q dynamo, in which the toroidal magnetic 
field is generated by the shear of the poloidal magnetic field by the 
differential rotation U,, which is also called the Q effect. The decrease 
in the rotation with radius and equation (1) imply that By should be 
anti-correlated with Bp. To complete the dynamo cycle, the poloidal 
magnetic field has to be generated by the toroidal electromotive force 
Ep whose main contribution is froma diagonal component agp, which 
is the so-called a effect. Eg is, therefore, correlated or anti-correlated 
to By depending on the sign of a,, (equation (2)). This complete cycle 
forms a dynamo | waye, that oscillates with a period of 
Prheory = 2 App 2k | (refs. 22,23) and propagates in the direc- 
tion app V O x ep according to the Parker-Yoshimura rule*“", Here k, is 
the wavenumber of the dynamo wave in the vertical direction. In this 
theoretical description, we have supposed that contributions from the 
other a; components and the turbulent resistivity tensor £; are sub- 
dominant. We will show that it is the case in the following. 

First, we confirm that the simulation set-up employed can capture 
the fastest-growing mode of the magnetorotational instability in the 
outer region of the hypermassive neutron star (Extended Data Fig. 2). 
Therefore, the turbulent state is developed. Figure 3 plots the butterfly 
diagrams for By, Bgand Egat R = 30 km. The top left and top right panels 
show the anti-correlation between By and Bp, which indicates the Q 
effect. To quantify the correlation between By and Ep, we compute 
the Pearson correlation C,(X, Y) between the two quantities X and Yin 
the bottom right panel of Fig. 3. This figure shows that By and &4 
are anti-correlated for z< 15 km, where the pressure scale height at 
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Fig. 4 | Electromagnetic signals from the remnants of the binary neutron 
star merger. Top left, Angular distribution of the Poynting flux ona sphere with 
r= 500 km at t- Cinerger “150 ms. Top right, Electron fraction distribution for 
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the magnetically driven post-merger ejecta (dashed) and the dynamical ejecta 
(dotted). Bottom, Same as the top right panel but for the terminal velocity (left) 
and entropy per baryon (right). 


R=30 kmis~14.6 km. The correlation between the electromotive force 
and mean current J; = (v x B) is small, and thus, the £; tensor can be 
L 


neglected (Extended Data Fig. 6). Therefore, the generation of the 
mean poloidal magnetic field is determined primarily by the a; tensor. 
Although a;r has a non-negligible contribution, the contribution of 
QB, dominates in the turbulent region for z < 10 km (Extended Data 
Fig. 6). This shows that a, is the main component, and it is plotted in 
the bottom right panel of Fig. 3. We also confirmed that the a20 dynamo 
is irrelevant compared to the «Q dynamo (Extended Data Fig. 7). 

With these quantities, we can predict the period of the «Q dynamo, 
as listed in Table 1. We measured the shear rate q = —d In Q/dInR at 
the selected radius and choose the wavenumber k, corresponding to 
the pressure scale height, which is the most extended vertical length 
in the turbulent region. The sixth and seventh columns report the 
predicted period of the «Q dynamo and the period measured in the 
butterfly diagram. The agreement at R = 30 kmis remarkable. We show 
the comparison at different radii from R = 20 to 50 kmin the table, and 
the agreement is also reasonable. In addition, since a,,is negative and 
qis positive, the dynamo wave propagates along the direction of the 
Parker-Yoshimura rule, that is, the z direction (http://www2.yukawa. 
kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/Br_core_cut_9km.m4v). 
With these findings, we conclude that the dynamo in our simulation 
can be interpreted as an «Q dynamo, and it builds up the large-scale 
magnetic field in the remnant of the hypermassive neutron star 
(Fig. 1). 

After the development of the coherent poloidal magnetic field 
due to the «Q dynamo and the resultant efficient magnetic winding, 
a magnetic-tower outflow is launched towards the polar direction 
(http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/ 
DD2_135_135 Dynamo.mp4; Methods). The mean poloidal magnetic 
flux is generated in the region where the magnetorotational instability 
is active. The mean magnetic field deep inside the core (p > 10" g cm°), 
which could be a relic of the initial field, stays buried belowr< 10 km 
in the polar region throughout the simulation (http://www2.yukawa. 
kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_Mean_Poloidal_ 
Flux.mp4; Methods). This indicates that the large-scale field generated 


by the «Q dynamois responsible for launching the jet. We confirm that 
the low-resolution simulation cannot capture the launching of the 
strong Poynting-flux-dominated outflow or the enormous post-merger 
ejecta” (see also Extended Data Fig. 8). 
The luminosity of the Poynting flux, which is defined by Lpoy = 

-$ soom V £ (Togn 42, where Tami is the stress energy tensor for 
the clectromagnetie field, is -10% erg s ‘at the end of the simulation 
when t- tmerger ¥ 150 ms. The high Poynting flux is confined in the region 
with 0 < 12° where @ is a polar angle, as plotted in the top left panel of 
Fig. 4 (see also Extended Data Fig. 8). The luminosity and jet opening 
angle are broadly compatible with some of the observed short, hard 
gamma-ray bursts”. We estimate the terminal Lorentz factor l. by the 
magnetization parameter 0,c = b*/4npc at the light cylinder radius 
Nic = ¢/2 x 40 km((2/7, 000 rad s-!) ‘contained inaregion@<12° 
(refs. 43,44). We found the baryon loading is still severe, for example, 
Mouttlow © 107? Mg s™ in the polar region. Nonetheless, .. increases 
with time but fluctuates. It reaches -10-20 at the end of the simulation. 
Therefore, if the magnetic reconnection efficiently dissipated the 
Poynting flux, a Lorentz factor of -10 could be possible“. Also, the 
evacuation due to the strong Poynting flux in the polar region continues 
at the end of the simulation. Therefore, a higher F. could be possible 
after the long-term evolution. 

The dynamical ejecta driven by the tidal force and the shock wave 
during the merger phase represent ~0.002 M, (ref. 34). After the devel- 
opment of the coherent magnetic field, we observe anew component in 
the ejecta that is driven by the Lorentz force, namely, the magnetically 
driven wind“. The mass of this component is -0.1 M, at the end of the 
simulation t — tnerger “ 150 ms. The electron fraction of the post-merger 
ejecta shows a peak around -0.2 with an extension to ~0.5. The terminal 
velocity of the post-merger ejecta peaks around ~0.08c, where cis the 
speed of light. Therefore, it could contribute to a bright kilonova emis- 
sion through the synthesis of r-process heavy elements”**®. 

We speculate that the high-luminosity state and the post-merger 
mass ejection would continue for O(1) s, which corresponds to the 
neutrino cooling timescale**“’”. As reported in refs. 27,34, the torus 
expands due to the angular momentum transport facilitated by the 
magnetorotational instability-driven turbulent viscosity after the 
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neutrino cooling becomes inefficient. As a result, the funnel region 
above the remnant neutron star expands and the jet could no longer be 
collimated. A long-term simulation for O(1) s ofa remnant of a massive 
neutron star is future work to be pursued. 

Also, a simulation with initially weakly magnetized binary neutron 
stars is necessary to confirm the picture reported in this article because 
the saturated magneto-turbulent state and the generation timescale 
of the coherent poloidal magnetic field due to the magnetorotational 
instability could depend on the initial magnetic field strength and 
flux*®. The recent large-eddy simulations of a magnetized binary neu- 
tron star merger may give a hint about this issue*””°. Even if the initial 
magnetic field has amuch weaker strength or different topology from 
the one we assumed, the saturation strength and field profile due to 
the Kelvin-Helmholtz instability in the merger remnant are like those 
that we found in this article. Also, how the mean poloidal magnetic field 
is set in reality after the merger is an open problem. If the mean poloidal 
magnetic field just after the merger is a relic of the pre-merger poloidal 
field, that is, 10-10" G at maximum, it may take In(10*) x 0.02s ~ 0.2s 
to reach the saturation strength of 10-10" G, as we assume the mean 
poloidal field is exponentially amplified with the period of the aQ 
dynamo” (Fig. 2 and Table 1). Consequently, the jet may be launched 
at O(0.1) s after the merger in reality. However, the interior structure 
of the magnetic field in the pre-merger neutron stars is not well under- 
stood, and the magnetic reconnection of the fluctuating poloidal field 
generated by the Kelvin-Helmholtz instability may enhance the mean 
poloidal field after the merger. This situation should be explored as a 
future task. 

To summarize, we generated a large-scale magnetic field 
in a long-lived remnant of a binary neutron star merger using a 
super-high-resolution neutrino-radiation magnetohydrodynamics sim- 
ulation in general relativity. The magnetorotational instability-driven 
aQ dynamo generates the large-scale magnetic field. As a result, the 
launch of the Poynting-flux-dominated relativistic outflow in the polar 
direction and an enormous amount of magnetically driven wind are 
induced. Our simulation suggests that the magnetar engine generates 
short, hard gamma-ray bursts and bright kilonovae emission, which 
could be observed in the near-future multi-messenger observations. 


Methods 

Numerical method 

Our code implements the Baumgarte-Shapiro-Shibata- 
Nakamura-puncture formulation to solve Einstein’s equation! *. 
The code also uses the Z4c prescription to suppress the constraint 
violation”. The fourth-order accurate finite difference is used as a 
discretization scheme in space and time. The sixth-order Kreiss—Oliger 
dissipation is also employed. The HLLD Riemann solver and the con- 
strained transport scheme” are used to solve the equations of motion 
of the relativistic magnetohydrodynamic fluid. The neutrino-radiation 
transfer is solved by the grey M1+GR-leakage scheme” to take into 
account neutrino cooling and heating. 


Initial data 

We adopted quasi-equilibrium initial data of irrotational binary neu- 
tron stars in the neutrino-free beta equilibrium derived in a previous 
paper” using the public spectral library LORENE (ref. 59). The ini- 
tial orbital angular velocity is Gm)Q)/c? = 0.028 with my = 2.7 M,, The 
residual orbital eccentricity is -10°. For our high-resolution study, the 
data are remapped onto the computational domain described in the 
‘Grid set-up’ section. 


Grid set-up. We employ the static mesh refinement with 2:1 refine- 
ment, that is, a coarser domain with a grid resolution twice that 
of a finer domain. All the domains are composed of concentric 
Cartesian domains with a fixed grid size N. The size of the grid is 
2N x 2Nx Ninthex, y and z directions, and we assume orbital plane 


symmetry. We employ 16 domains with N = 361 and the finest grid 
resolution of AX finest = 12.5 m. The first three finest domains, whose 
sizes are [-4.5: 4.5 km]? x [0: 4.5 km], [-9: 9 km]? x [0: 9 km] and 
[-18: 18 km]? x [0:18 km], respectively, are employed to resolve the 
Kelvin-Helmholtz instability, which emerges on a contact interface 
(shear layer) when the two neutron stars merge. The initial binary 
separation is ~44 km, and the coordinate radius of the neutron star is 
10.9 km. Thus, the fourth finest domain with [-36: 36 km]? x [0:36 km] 
covers the entire binary neutron star. 

We started each simulation with AXnes = 12.5 mand 16 static mesh 
refinement domains. At ~30 ms after the merger, we removed the two 
finest domains with AX finest = 12.5 m and Ax = 25 m and continued the 
simulation with AX finest = 50 m until -50 ms after the merger. Then, we 
removed the domain with AX finest = 50 mand continued the simulation 
with AXinest = 100 m. The timing for removing the finer domains was 
determined by monitoring the magnetorotational instability quality 
factor so as not to go below the critical value of 10 after the removal (see 
below). This strategy allowed us to capture the efficient amplification 
of the magnetic field through the Kelvin-Helmholtz instability and 
resolve the magnetorotational instability inside the remnant while 
saving computational costs. 

To check the validity of removing finer domains, we continued 
the simulation with AXfines: = 12.5 M UP to t — tmerger * 33 ms (See the blue 
dashed curve in Fig. 1). We observed an -10% decrease in the poloidal 
electromagnetic energy in the simulation with AXines, = 50 m com- 
pared to that with AXfnest = 12.5 m around this time. Nonetheless, the 
poloidal electromagnetic energy increased around t- t merger * 40 ms 
in the simulation with AX finest = 50 m (the green dashed curve in Fig. 
1) due to the a effect we discuss in the main text. Also, there was no 
substantial decrease in the toroidal electromagnetic energy after 
removing the finer domains at t- tmerger “ 30 ms. Similarly, there was 
no visible deterioration after the second removal of the finer domain 
at t- t merger “ 50 ms. 

For the convergence test, we also performed the simulations with 
AXines = 100 (200) m and N = 361 (185). The number of domains was 13. 


Equation of state. We extended the original DD2 equation of state to 
the low-density and temperature region with the Helmholtz equation 
of state®’. Because any high-resolution shock-capturing scheme does 
not allow a vacuum state, we employed the atmospheric prescription 
outside the neutron stars. Specifically, we set a constant atmospheric 
density of 10? g cm” inside r < L,,,, and assumed a power-law profile 
Pam = Max[103(Ljtm/N) gem, Pq] for r> Lam Where Lam = 36 km. The 
Helmholtz equation of state employed determines the floor value of 
the rest-mass density, which is pa = 0.167 g cm”. We also assumed a 
constant atmospheric temperature of 10° MeV. 


Kelvin-Helmholtz instability 

The top panel of Extended Data Fig. 1 is the same as the inset in Fig. 1 
but with additional data. The solid blue, solid brown and solid purple 
curves are the results with AX ines = 12.5, 100 and 200 m, respectively, 
with Bo max = 105° G. The solid red and solid orange curves show the 
results with AX finest = 12.5 m with By max = 10% and 10" G, respectively. 
The cyan curve shows the result with AX finest = 18.75 mand Bo, max = 10“G. 
The dotted blue and dotted red curves plot the result with 
AX finest = 12.5 m and Bo max = 10 G magnified by the square of the ratio 
of the initial magnetic field strength, that is, (10155/1017 =10° and 
(105/10!4)" = 107, respectively. 

This figure suggests the following two points. First, the solid and 
dotted blue curves and the solid and dotted red curves overlap until the 
back reaction starts to activate, which happens when the electromagnetic 
energy reaches -3 x 10” erg. This implies that the magnetic field is passive 
for t- Cmerger £ 1 ms, which is the linear phase. The saturation of the elec- 
tromagnetic energy due to the Kelvin-Helmholtz instability is likely to 
be -1-3 x 10% erg corresponding to O(0.1)% of the internal energy?°?*!, 
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The second point is that the growth rate of the electromagnetic 
energy in the linear phase depends on the grid resolution not on the 
initial magnetic field strength. To quantify the dependence of the 
growth rate on the grid resolution and initial magnetic field strength, 
we estimated the growth rate by fitting the electromagnetic energy as 
Emag(t) = A exp (Oxu(t — tmerger)) for O $ t- tmerger £ 1 ms, which corre- 
sponds to the linear phase. The bottom panel of Extended Data Fig. 1 
plots the estimated growth rate as a function of the initial magnetic 
field strength. The symbols denote the grid resolution. With 
AXinest = 12.5 m, the growth rate was -7 ms“, irrespective of the initial 
magnetic field strength. The growth rate increased with the grid resolu- 
tion. With Axfnes: = 100 or 200 m, the efficient amplification cannot be 
captured. All the properties are consistent with those of the Kelvin- 
Helmholtz instability, that is, the smaller the scale of the vortices is, 
the larger the growth rate is. The reality is located in the leftward and 
upward direction in this diagram, that is, the weak magnetic field 
observed in binary pulsars (the leftward direction in the diagram) and 
infinitesimal grid resolution (the upward direction in the diagram). 
Therefore, we anticipate that the electromagnetic energy will saturate 
in avery short timescale, «1 ms after the binary neutron star merger, 
irrespective of the magnetic field strength at the onset of the merger 
unless a black hole is promptly formed®*>*, 

We estimated how small AXjines: Should be to simulate the amplifica- 
tion due to the Kelvin-Helmholtz instability starting from the ‘realistic’ 
initial magnetic field strength of 10" G (ref. 36). The saturation strength 
found inthe simulation with AX finest = 12.5 m is likely to be |B] ~ 10165 G. 
We fitted the growth rate as a function of the grid resolution from 
Extended Data Fig. 1. It was Oęķu(ms™”) = 90/AXines(M), where we assumed 
the inverse proportionality would originate due to the Kelvin-Helm- 
holtz instability”. Assuming that the initial magnetic field strength 
was 10” Gand that the lifetime of the shear layer was 2 ms, we estimated 
the required growth rate as 


1 ( 1016-5 G 


2 
= Z -1 
TT 101G ) 12.7 ms™!. (3) 


Oxy 


Therefore, the required grid resolution is 


90 
AXfinest = Ban ~ 7.1m. (4) 


Magnetorotational instability, neutrino viscosity and 
neutrino drag 

To quantify how well the magnetorotational instability is resolved 
in our simulation, we estimated the rest-mass-density-conditioned 
magnetorotational instability quality factor, as defined by 


Amri)p 1 S,Ame dx 
Ax Ax fax ’ 


(Quri)p = (5) 


where Ayrı = 2087/ (Y4rp 2) is the fastest-growing mode of an ideal 
magnetorotational instability. As the condition in terms of the 
rest-mass density, we define a remnant core and a remnant torus as 
fluid elements with p > 10" and <10”? g cm”, respectively. Furthermore, 
we exclude the core region above 10** g cm” in the estimate of the 
quality factor because in sucha region, the radial gradient of the angu- 
lar velocity is positive, and it is not subject to magnetorotational insta- 
bility”, as plotted in the top panel of Extended Data Fig. 2. Also, we 
introduce a cutoff density of 10’ g cm” for the torus to suppress the 
overestimation of the quality factor in the low-density polar region. 
The middle and bottom panels of the figure show that the 
fastest-growing mode of the magnetorotational instability is well 
resolved, both in the core and torus, throughout the simulation. Con- 
sequently, magneto-turbulence develops and is sustained. However, 
the magnetorotational instability, particularly inthe core throughout 


the entire simulation and in the torus for t — tnerger £ 80-90 ms, cannot 
be resolved when AXfinest = 200 m. Thus, the turbulence is not fully 
produced in sucha low-resolution run. 

One caveat is that the neutrino viscosity and drag could affect the 
magnetorotational instability as diffusive and damping processes”. The 
former (latter) becomes relevant when the neutrino mean free path is 
shorter (longer) than the wavelength of the magnetorotational instabil- 
ity. Given a profile of the merger remnant, suchas the density p, angular 
velocity Q, temperature Tand hypothetical magnetic field strength 8z y” 
we solve the two branches of the dispersion relations to quantify the 
neutrino viscosity and drag effects“: For the neutrino viscosity: 


(0+ Enor E] +e? +@]-48 =0, 6) 


and for the neutrino drag: 


[(o+ ore] +Ë [7 + È| -4Ř =0, (7) 


where 6 = 0/2, k = kval 2, R = P/2,¥, = WO/X and Ñ, = T,/2. 0 and 
k are the growth rate and wavenumber of the unstable mode of the 
magnetorotational instability. vą = Bi, ,/V4npis the Alfvén wave speed. 
r is the epicyclic frequency. v, and F, are the neutrino viscosity and 
drag damping rate, respectively. Reference © provides fitting formulae 
for them as functions of the rest-mass density and temperature: 


2 
vy =12x 10°( ) cm? s7}, (8) 


Pp ( T 
108 g cm7? 10 MeV 


(9) 


T \6 
= Bf \ s 
I,=6x10 (om) sv. 


Also, the neutrino mean free path J, is fitted by 


=f 2 
T 
L=2x10| —? ( ) l 
y x (= g cm~? ) 10 MeV cm 


Extended Data Fig. 3 plots the growth rate of the magnetorota- 
tional instability for a given profile of a remnant of a massive neutron 
star and a hypothetical value of Bisp We use the profiles on the orbital 
plane at t- tmerger “ 10 ms. The purple curve denotes the boundary where 
the condition ¥, or f, x 1 is met®. Inside the boundary, the neutrino 
viscosity or the neutrino drag substantially suppresses the growth rate. 
Outside it, the growth rate is essentially the same as for the ideal mag- 
netorotational instability. Above it, we plot the azimuthally averaged 
magnetic field strength 54m obtained in the simulation. Because of the 
efficient amplification due to the Kelvin-Helmholtz instability just 
after the merger, the neutrino viscosity and drag effects are irrelevant 
inthe entire region of the merger remnant. 


(10) 


Generation of mean poloidal magnetic flux due to the 
magnetorotational instability 

To validate our interpretation that the «Q dynamo is responsible for 
generating the large-scale magnetic field, we evaluate the mean poloi- 
dal magnetic flux on a certain sphere of radius r: 


m/2 


Pg (N= 2f B,r? sin 0 dé, (11) 
0 


where 8B, denotes the radial mean field in the spherical-polar 
coordinate. 

Extended Data Fig. 4 plots the evolution of the mean poloidal 
magnetic flux. We selected 8 and 20 kmas representative radii for the 
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magnetorotational instability inactive and active regions (see the top 
panel of Extended Data Fig. 2). Onthe one hand, inthe simulation with 
AX finest = 12.5 m, the mean poloidal magnetic flux on the sphere with 
r=8 kmexhibits intensive time variability for t — tmerger £ 20 ms, which 
reflects the oscillations of the remnant core, that is, compression and 
decompression. The poloidal flux gradually decreases, presumably 
due to the magnetic reconnection, whichis imprinted as the decrease 
of the total poloidal field energy for 10 $ t- tmerger $ 20 ms in Fig. 1. 
Nonetheless, it relaxes to a roughly constant value and stays there until 
t-tnerger S 120 ms. This reflects that magnetorotational instability is 
inactive in this region, that is, the mean poloidal magnetic flux is not 
generated. On the other hand, it is evident that the mean poloidal flux 
on the sphere with r= 20 km is generated around t- tmerger ¥ 25 ms in 
the high-resolution simulation. Such efficient generation of the mean 
poloidal flux is not observed in the simulation with AXx¢,es: = 200 m for 
t- tmerger £ 80 ms, as plotted with the green dashed curve in the figure 
(but see below for details about the slight generation of the mean 
poloidal flux for 40 $ t — tmerger £ 60 MS). FOr t — t merger z 100 ms, the 
mean poloidal magnetic flux is generated even in the low-resolution 
simulation because the magnetorotational instability starts to be 
partially, not fully, resolved in the low-density region (see the middle 
and bottom panels of Extended Data Fig. 2). Therefore, we conclude 
that the magnetorotational instability is responsible for generating 
the mean poloidal magnetic flux. 

Extended Data Fig. 5 plots meridional profiles of the mean radial 
magnetic field B, at selected time slices for the simulation with 
AX finest = 12.5 (200) m in the left (right) column (see also the visualiza- 
tion: http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/ 
SAKURA/movie_Mean_Poloidal_Flux.mp4). Itis evident that the mean 
poloidal field is generated in the region where the magnetorotational 
instability is active (p < 10° g cm”), particularly, in the low-latitude 
region 0.24 < 0 < 0.5m. Then, it propagates towards the high-latitude 
region and generates the jet. We also point out that at a radius of 
9-10 km, the polar region has a weak magnetic field. This indicates that 
the poloidal field below this region does not contribute towards the 
launch of ajet and stays buried throughout the simulation. 

The low-resolution simulation also shows some amplification of 
the mean poloidal field in the low-density region with p < 10” g cm”, 
where the magnetorotational instability is resolved. (see the visuali- 
zations for the mean poloidal flux and magnetorotational instability 
quality factor in the simulation with AX nese = 200 m: http://www2. 
yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_Mean_ 
Poloidal_Flux_Low.mp4 and http://www2.yukawa.kyoto-u.ac.jp/-kenta. 
kiuchi/anime/SAKURA/movie_Mean_MRI_Qfac_Low.mp4). It shows 
qualitatively similar behaviour to the simulation with AX finest = 12.5 m, 
but the efficiency for generating the mean poloidal magnetic field is 
much lower, which leads to a weaker Poynting flux luminosity (see 
the ‘Detailed properties of the Poynting-flux-dominated outflow and 
magnetically driven post-merger ejecta’ section). 


aN dynamo 

To understand the dynamo process inthe simulation, we used the mean 
field theory with an axisymmetric average. We consider an axisym- 
metric average because it corresponds to the symmetry of the differ- 
ential rotation in the system. In the mean field dynamo theory, we 
assume that the velocity field U and the magnetic field B are, respec- 
tively, composed of the mean velocity field Uand velocity fluctuations 
usothat U = U + uand of the mean magnetic field B and the magnetic 
fluctuations b. We then average the induction equation, which gives 
equation (1), where ē = ux b is the electromotive force due to the 
fluctuations. To close the system, the electromotive force is often 
expressed as a function of the mean magnetic field, and its spatial 
derivatives as described by equation (2). a, and £; are, respectively, 
tensors expressing the contributions of the mean magnetic field B and 
its derivatives J = V x B to the electromotive force. These tensors 


should not depend on B. First, we focus on how the mean poloidal field 
is generated and then look at the toroidal field. 

To show how the mean poloidal field is generated by the alpha or 
beta effect, we assume that the mean velocity field is composed of the 
rotation speed U = Re, and project the averaged induction equation 
(1) in the radial direction in cylindrical coordinates, which gives: 


ðBr = -d,|(U xB), + Ep] = 0.89 = -0, (a19)B; + By(V xB),). 2) 


Similarly, by projecting the averaged induction equation inthe vertical 
direction, the generation of the mean vertical field B, is due to the 
azimuthal component of the electromotive force €,. For the generation 
of the mean toroidal field, the Q effect (that is the winding of the mag- 
netic field by differential rotation) is also important and must be com- 
pared tothe contributions of both the radial and vertical components 
of the electromotive force &, and &,. To estimate which of the mean 
magnetic field 8; or the derivatives of the mean magnetic field J, con- 
tributes the most to the electromotive force €;, we compute the Pearson 
correlation coefficients between the quantities €,and Y, = B,or J, with 
the following formula: 


S Ei — EG — V) dt 


Cp(E, Y;) = 
VUE- Ene S- HY do) 


i (13) 


where (-), represents a time average. 

Inthe following sections, we present acomplementary analysis of 
the other contributions to the generation of the mean magnetic field 
than the apg effect and Q effect. We, therefore, show the other correla- 
tions between the electromotive force and the magnetic field and the 
estimated values of the tensor component. Several methods can be 
used for the values of the alpha tensor components. The simplest one 
is to estimate them from the correlations”, but this method assumes 
that onecomponentis dominant. To take into account the off-diagonal 
contributions, we compute the values of the alpha tensor coefficients 
in this study by using a singular value decomposition ina least-squares 
fit of the mean current data and mean field’. 


Generation of the mean poloidal field Bz. As shown in equation (12), 
the generation of the axisymmetric poloidal field is due to the curl of 
the toroidal component of the electromotive force in the averaged 
induction equation. In this section, we show the correlations between 
the toroidal electromotive force €, and the magnetic field components 
B,and the mean current J; (top and middle panels of Extended Data Fig. 
6). The low correlations with the mean current J, shows that its contri- 
butions to €, (that is the £; tensor) can be neglected. In the same way, 
the vertical magnetic field contribution 8, canbe neglected. Extended 
Data Fig. 6 also shows the anti-correlation of 8, and ë, and also that 
the radial magnetic field Bgis correlated to &y. This is because the radial 
field Bg is anti-correlated to the toroidal magnetic field By due to the 
Q effect (Fig. 3). Since the mean toroidal field B, is anti-correlated to 
the toroidal electromotive force €,, the mean radial field Bẹ is corre- 
lated to it. To confirm that the contribution of B, dominates, we first 
compute the non-diagonal components of the alpha tensor (bottom 
panel of Extended Data Fig. 6). We then compare the contribution of 
B, and Bp, that is, respectively, the time-averaged values of a,,B, and 
pB, in the first 10 km: 


(agpBo), 


—— = 1.87. 
(AprBr), 


(14) 


Generation of the mean toroidal field. In the main text, we show that 
the Qeffect is important for the generation of the mean toroidal field. 
In this section, we check whether the contribution of the a effect 
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through the poloidal components of the electromotive force Eg and &, 
is substantial, in which case the dynamois called an @°Q dynamo instead 
of an «Q dynamo. First, we confirm that the correlations between the 
electromotive force and the mean current J, (right panels of Extended 
Data Fig. 7) are lower than with the mean magnetic field 8; (left panels 
of Extended Data Fig. 7). The contributions from the mean current can, 
therefore, be neglected. For the mean magnetic field, some compo- 
nents, for example, the mean toroidal field By, are strongly correlated 
with the radial component of the electromotive force &p. This indicates 
that the a effect might be important to the generation of the mean 
toroidal field. To compare the strength of these two effects, the a effect 
andthe Q effect, we computed the corresponding a tensor components 
ær; and a,; and estimated the ratio of the two dynamo numbers 
Cx = Max(|agil, |a,;))R/n for ic [R, 6, z] and Cy = OR?/n, where 7 is the 
resistivity, in the turbulent region averaged for one scale height, which 
gives atR=30km: 


GA OR 
Cy  mMAaX(|aril, lazil) 


30.8. (15) 


The contribution of the a effect to the generation of the mean toroidal 
field can, therefore, be reasonably neglected as the Q effect dominates 
the generation of the mean toroidal field. The dynamo in our simulation 
can, therefore, be interpreted as an «Q dynamo. 


Detailed properties of the Poynting-flux-dominated outflow 
and magnetically driven post-merger ejecta 
The link http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/ 
SAKURA/DD2_135_135 Dynamo.mp4 isa visualization for the rest-mass 
density (top left), the magnetic field strength (top, second from left), 
the magnetization parameter (top, second from right), the unbound- 
edness defined by the Bernoulli criterion (top right), the electron 
fraction (bottom left), the temperature (bottom, second from left), 
the entropy per baryon (bottom, second from right) and the geodesic 
criterion (bottom right) ona plane perpendicular to the orbital plane. 
The top panel of Extended Data Fig. 8 plots the angular distribution 
of the luminosity of the Poynting flux, which is calculated as follows: 


ayer? (T gem) A, (16) 


Lpoynting() = -f 
rx500km 

where wand Ware the lapse function and the conformal factor, respec- 

tively. The high luminosity of -2-8 x 10” erg s” angle ‘is confinedina 

region with 0 < 12°. 

The middle panel plots the jet-opening-angle-corrected luminos- 
ity and luminosity of the Poynting flux (green and blue, respectively) 
as functions of the post-merger time. According to refs. 67-69, the 
luminosity of the Poynting-flux-dominated outflow, which is driven by 
the efficient magnetic winding from the remnant of the binary neutron 
star merger, can be estimated as 


Lpoy = 
(17) 


S Lpoynting(9) sin 6d6 = 10°% erg si mn 1 Ta j (i s-1 ) ý 


where 8, is the azimuthally averaged poloidal magnetic field. In the 
current simulation, itis -10” G. Thus, the Poynting flux luminosity found 
in the simulation is consistent with this estimation. The 
jet-opening-angle-corrected luminosity is ~10” erg s”. Thus, ifwe assume 
that the conversion efficiency to gamma-ray photons is ~10%, it is com- 
patible with the observed luminosity of short, hard gamma-ray bursts”. 

The bottom panel plots the ejecta as a function of the post-merger 
time. The solid curve denotes the result from using the Bernoulli cri- 
teria. The coloured region in the inset shows the error of the baryon 
mass conservation, whichis below 10” M, throughout the simulation. 


Convergence study and effect of the initial large-scale 
magnetic field 
Because we chose an initial large-scale magnetic field strength that is 
muchstronger than those observed in binary pulsars*, we need to clarify 
whether suchastrong large-scale field is responsible for launching the jet. 
We begin by estimating the magnetic winding timescale origi- 
nating from the initial magnetic field. Suppose we consider a binary 
neutron star merger witha highly magnetized end of 10" G, as observed 
for pulsars. In that case, the timescale for the magnetic winding to 
reach saturation is 


By x 10165 o( (18) 


Br )( Q ($) 
5x10" G/\8,000rads-1/\8s/’ 


where we assume that the compression due to the merger amplifies 
the initial poloidal field by a factor of five, which should be propor- 
tional to p°” because of conservation of the magnetic flux. We also 
assume a Saturation field strength of 10° G, as suggested by the 
super-high-resolution simulation (Extended Data Fig. 1). Therefore, 
the magnetic winding originating from the initial magnetic field should 
be minor or irrelevant in reality. 

However, the magnetic winding timescale is much shorter 
than those in reality in both simulations with AX fines, = 12.5 or 200 m 
because of the assumed strong initial field. The low-resolution mesh 
is fine enough to resolve the compression and winding but can only 
partially capture the Kelvin-Helmholtz and the magnetorotational 
instability. Therefore, there would be no striking difference between 
the simulations for the jet launching mechanism if such a strong ini- 
tial poloidal field and subsequent magnetic winding were respon- 
sible for it. The middle panel of Extended Data Fig. 8 shows that the 
Poynting-flux-dominated outflow is launched at t — tnerger * 60 ms 
and reaches -10° erg s7 in the low-resolution simulation with 
AX finest = 200 m. However, in the same plot, the super-high-resolution 
simulation with AX finest = 12.5 m shows that the Poynting-flux-dominated 
outflow is launched at t- tmerger “ 35 MS. The luminosity reaches 
~10™ erg s”. The difference is striking. As discussed in Extended Data 
Fig. 5, the efficiency of the generation of the mean poloidal flux in 
the super-high-resolution simulation is much higher than that in the 
low-resolution simulation. This indicates that the efficient «Q dynamo 
is responsible for launching the strong jet. 

The bottom panel of Extended Data Fig. 8 also suggests that the 
Lorentz force-driven post-merger ejecta is launched in the run with 
AX finest = 200 m. However, the ejecta mass is about one order of mag- 
nitude smaller than in the run with AXinest = 12.5 m, showing that the 
enhanced activity of magnetohydrodynamics effects by the efficient 
dynamo action plays an important role in ejecting matter”. 


Data availability 

The raw simulation data (265TB) that support the findings of this study are 
available from the corresponding author upon request. The prerequisite 
for a data request is having a server where the data can be transferred to. 


Code availability 
The numerical relativity code and data analysis tool are available from 
the corresponding author upon reasonable request. 
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Extended Data Fig. 1| Kelvin-Helmholtz instability growth rate at the merger. respectively. (Bottom) Dependence of the growth rate of the electromagnetic 


15.50 


(Top) Same as the inset in Fig. 1 in the main article, but with AX fines: = 18.75 m for energy at the merger due to the Kelvin-Helmholtz instability on the initial 
Bo,max = 10“ G (cyan). The blue- and red-dotted curves show the evolution for magnetic field strength Bo max and grid resolution. The error is due to the fitting. 
AXpinest = 12.5 mand Bo max = 10! G magnified by a factor of 10° and 10°, Data are presented as mean median values +/- SD. 
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Extended Data Fig. 2 | Magnetorotational instability inside the merger 
remnant. (Top) The radial profile of the angular velocity (blue) and the rest-mass 
density (green) on the equatorial plane at t- frerger = 50 ms. Magnetorotational 
instability is inactive in a region with x < 9 km and p z 105 g/cm*. The inset shows 
the same profiles with the logarithmic scale in the radial direction. (Middle) 
Magnetorotational instability quality factor in a core region as a function of 
time. The remnant core is defined by a region with p210” g/cm’. The blue curve 
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denotes the employed finest grid resolution of 12.5 m. At t - t merger * 30 ms, the 
two finest domains with AXfnest = 12.5 m and AXfnest = 25 m are removed. Thus, 

the employed grid resolution is AX fines: = 50 m plotted with the green curve. 

Att- tmerger “ 50 ms, the finest domain with AXxgines: = 50 m is removed and the 
subsequent evolution with AXfines: = 100 mis plotted with the cyan curve. The 
result with AX fines: = 200 m is plotted with the purple curve. (Bottom) The same as 
the middle panel, but for a torus defined by a region with 10’ g/cm?<ps10" g/cm. 
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Extended Data Fig. 3 | Neutrino viscosity and drag on the magnetorotational growthrate is essentially the same as the ideal magnetorotational instability. The 
instability. Growth rate of the fastest growing mode of the neutrino viscous/drag red curve denotes the z-component of the azimuthally averaged magnetic field 
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magnetorotational instability as a function of the radius of the remnant massive strength in the simulation. The orange-solid (dashed) curve denotes the 
neutron star and the hypothetical value of the z-component of the magnetic field three-dimensional data for the z-component on the x = 0 axis in the simulation 
Biyp We take the simulation data on the orbital plane at t — tmerger = 10 ms. The With AX finest = 12.5(200) m. 
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Extended Data Fig. 4| Generation of the mean poloidal-magnetic flux inside r=8and20km, respectively, which are representative for the magnetorotational 
the merger remnant. Mean poloidal-magnetic flux as a function of the time instability inactive and active region, in the simulation with AX fines: = 12.5 m. The 
after the merger. Blue and green-solid curves denote the flux on the sphere of dashed counterparts denote the simulation with AX ines, = 200 m. 
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Extended Data Fig. 5 | Generation of the mean radial magnetic field. 
Meridional profile of the mean radial magnetic field, B,, at ~ 10 ms (top), 
30 ms (center), and 50 ms (bottom) for the simulation with AX finest = 12.5 m 
(left column) and 200 m (right column). The rest-mass density contour 
curves are also plotted. The visualization is the following link 


(http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_ 
Mean_Poloidal_Flux.mp4) for the super-high resolution simulation, and 
http://www2.yukawa.kyoto-u.ac.jp/-kenta.kiuchi/anime/SAKURA/movie_ 
Mean_Poloidal_Flux_Low.mp4 for the low-resolution simulation. 
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Extended Data Fig. 6 | a effect vs turbulent resistivity in the aQ dynamo. (Top) correlations between the toroidal electromotive force €, and the mean current 
Time-averaged correlations between the toroidal electromotive force Eg and the 


jj components for R = 30 km. (Bottom) Alpha tensor components ag; with 
mean magnetic field components B; for R = 30 km. (Middle) Time-averaged ie[R,z] forR=30km. 
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Extended Data Fig. 7 | aQ vs a’O dynamo. (Left) Time-averaged correlations between the poloidal electromotive force Ep (top) and £, (bottom) and mean magnetic 


field components 8; for R = 30 km. (Right) Time-averaged correlations between the poloidal electromotive force Ep (top) and €, (bottom) and mean current 
components J; for R=30km. 
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Extended Data Fig. 8 | Properties of the Poynting flux-dominated outflow plots the luminosity for the simulation with AX nes: = 200 m. (Bottom) Ejecta as a 
and post-merger ejecta. (Top) Angular distribution of the luminosity of function of the post-merger time. The solid curve denotes the ejecta satisfying 
the Poynting flux at the end of the simulation of t - t merger = 150 ms. (Middle) the Bernoulli criterion. The colored region in the inset shows the violation of 
Luminosity for the Poynting flux as a function of the post-merger time. The green the baryon mass conservation. The blue-dashed curve plots the ejecta for the 
curve is the jet-opening-angle corrected luminosity. The blue-dashed curve simulation with AX finest = 200 m. 
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